The research on exotic phases and the phase transitions has always been one of the central topics in statistical and condensed matter physics.^[1] A typical example is the continuous XY model on the square lattice, which has attracted much attention since Kosterlitz and Thouless (KT),^[2,3] discovered a phase transition at finite temperature without any symmetry breaking, obviously beyond Landau’s symmetry breaking theory for continuous phase transitions. Its low temperature phase, characterized by the bound vortex–antivortex pairs, is called the quasi-ordered or topological phase, with a divergent correlation length. The transition to the high temperature disordered (or paramagnetic) phase is attributed to the unbinding of these topological pairs.^[2,3]

Its discrete version, the q-state clock model exhibits even richer features in phase transitions, dependent on q, and arouses extensive interests and debates consequently. As is well known, if q ⩽ 4, the model has only one order-disorder transition of Ising type; when q approaches infinity, it then transits to the continuous XY model, with one KT transition.^[4] Furthermore, it is agreed that there exists a critical q_c, once q ⩾ q_c, the system has two transitions, separated by the topological KT phase in-between. Apparently, the discrete symmetry plays an important role on the transitions. There has been a large amount of studies on this model with different q’s.^[5–18]

So far, still there are debates on the very value of q_c, though many believed to be 5. Besides, the precise determination of the transition temperatures and types is also inconclusive. An early renormalization group (RG) analysis suggested two KT-type transitions^[19] for q = 5 case. Another recent density matrix renormalization group (DMRG) calculation^[20] with maximum size L = 256, estimated these two KT transition temperatures as T_c1 = 0.914(12) and T_c2 = 0.945(17) respectively. While, Monte Carlo (MC) simulations even gave inconsistent conclusions among themselves,^{[10,13,15,21,22]} mainly due to a different usage of the boundary twist and the physical quantities investigated to characterize the transitions, details in Refs. [21], [22], and [13]. For the transition temperatures, Borisenko et al.^[15] obtained T_c1 = 0.90514(9) and T_c2 = 0.95147(9) with size up to L = 1024, by calculating the susceptibility; and Kumano et al.^[21] predicted T_c1 = 0.908 and T_c2 = 0.944 using largest size L = 256, with investigation of the helicity modulus. All methods deal with relatively small system sizes.

In recent years, the tensor network states method has developed rapidly and become one of the most powerful numerical tools to study the phase transitions in both the classical and the quantum many-body systems.^[23–29] The tensor renormalization group method based on the higher-order singular value decomposition (abbreviated as HOTRG),^[28] has been applied successfully to the continuous XY model,^[30] the 3-dimensional Ising model^[28] and the Potts model.^[31,32]

In this article, we have investigated the phase transitions of the ferromagnetic five-state clock model on the square lattice by the HOTRG method. Firstly, by inspecting the specific heat and the spin–spin correlation, we confirmed that this model indeed has two phase transitions, separated by a narrow interval characterized with a power law correlation, apparently different to the low- or high-temperature Ising-type phases. Then we investigated the magnetic susceptibility, with performing a bond dimension scaling, only the upper transition temperature is obtained, at T_c2 = 0.9565(7). Finally, by studying the property of the fixed-point tensor from the RG flow, these two transition temperatures are located precisely at T_c1 = 0.9029(1) and T_c2 = 0.9520(1) respectively. The results agree quite well with the MC^[15,21] and the DMRG^[20] predictions.

The article is organized as follows: Section 2 introduces the two-dimensional q-state clock model and the numerical method we used. The results and discussions are presented in Section 3. At last, a brief summary is given in Section 4.

The Hamiltonian of the ferromagnetic q-state clock model with an external magnetic field is defined as 1 where ⟨ij⟩ denotes the summation over the nearest neighbors and θ_i is the angle of the spin at site i, which is one of the discrete set θ = 2πk/q, k = 0, 1, 2, …, q − 1. J is the coupling constant between neighboring spins, which is set to 1 for convenience. h is the applied weak magnetic field in unit of J/μ, where μ is the magnetic moment of each spin also set as 1.

For a classical statistical system, the tensor network states method usually starts by expressing the partition function as a product of all local tensors 2 where indices (l_i, r_i, u_i, d_i) denote the four legs (left, right, up, down) of the tensor at site i. And each local tensor T is defined by 3 where W’s come from the singular value decomposition of the Boltzmann factor, 4 As sketched in Fig. 1 of Ref. [28], the HOTRG coarse-graining scheme is a real-space renormalization, with each step reducing the system size by half. Iterating the process along the horizontal and the vertical directions alternatively, we can finally obtain the partition function and the free energy as well.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Temperature dependence of the internal energy and the specific heat with h = 0 and D = 40. Two obvious peaks indicate two possible phase transition from the specific heat.

We should again note that the HOTRG is credited for handling big size systems. Only log₂N steps are needed for the full contraction, given size N, which then could be very large to approximate the thermodynamic limit, and thus avoid the error inherent in the finite size scaling. In this work, we run the coarse-graining process until each physical quantity has converged. Normally, it takes 40 to 60 steps, so the system size is 2⁴⁰ to 2⁶⁰, approaching the infinity. One more thing should be mentioned is the size of the local tensor’s legs, or the bond dimension denoted as D, is equal to q initially. During each RG step, it is expanded quickly, then truncated to guarantee the process manageable. Usually, bigger the bond dimension, better the accuracy.

As illustrated in Fig. 1, we first calculate the internal energy and the specific heat. In the specific heat, there are two obvious peaks around 0.8 and 1.0, indicating two possible transitions. To be more explicit, we compute the spin-spin correlation function along the horizontal direction at three representative temperatures, which is expressed as 5 As shown in Fig. 2, their behavior is clearly different, indicating two phase transitions separating the ordered and the disordered Ising-type phases, as also suggested in Refs. [4], [6], [7], and [19].

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Spin–spin correlation along the horizontal direction with D = 40. Three different behaviors, each characterizing one phase respectively, can be distinguished clearly.

The specific heat cannot be used to determine the transition points, for which is derivatively continuous, as in the classical XY model on the square lattice, and thus whose peaks do not correspond to the transition temperatures. We then turn to investigate the system’s response to a static weak external magnetic field. We first evaluate the magnetization per spin, as shown in Fig. 3. Because of the discrete orientations of the spin, the usually used magnetization exhibits significant fluctuations. Thus, the definition of the magnetization used in Refs. [14], [15], [11], [12], [10], and [13] is adopted 6 Figure 3 shows the temperature dependence of the magnetization with two tentative magnetic fields. In the low temperature, the system is in the ordered phase with the magnetization approaches to 1. In the high temperature, the system is in the disordered phase with the magnetization equals to 0. In the critical phase, the magnetization takes values between 0 and 1. This behavior is similar to XY case.^[30]

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Magnetization with existence of weak external magnetic fields, and the corresponding magnetic susceptibility as defined in Eq. (7) with D = 40.

From the magnetization, we can calculate the magnetic susceptibility 7 to characterize the phase transition, as demonstrated also in Fig. 3, where one can observe a clear peak, signifying its analogy to the KT transition in the XY model. The peak appears around 0.999 with the applied field about 6 × 10⁻⁵. The magnetic susceptibility decays to zero exponentially above the peak temperature, indicating a phase transition occurs. However, the lower one does not manifest itself in the same way, maybe due to other transition mechanisms instead of KT type, which requires further investigations.

Following the procedure in Ref. [30], for a given D, we can determine the transition temperature by locating the peak positions of the magnetic susceptibility with respect to the external fields, and then extrapolating to the zero-field limit. As only one peak shows up, the upper transition between the quasi-ordered and the disordered phases is then predicted at T_c2 = 0.9561(10) by a power law fitting for D = 40, as presented in Fig. 4, 8 where a = 0.5666(310) and b = 0.2677(84). Here, we can see the behavior is very similar to the KT transition in the XY model,^[30] indicating it is also a KT-type transition. Furthermore, we perform the same process with respect to different bond dimensions to get the converged upper transition temperature, i.e., T_c2 = 0.9565(7) as shown in Fig. 5.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) The peak position of the magnetic susceptibility with respect to the external field, and a power law fitting is also performed for extrapolation to zero-field limit. The critical point is obtained as Tc2 = 0.9561(10) for D = 40 case.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) The upper phase transition temperature with respect to the bond dimension D by the magnetic susceptibility calculation. The converged transition temperature is Tc2 = 0.9565(7).

For a complete description of the transition temperatures, we turn to focus on the local fixed-point tensor. At each coarse-graining step, an optimized isometric matrix is obtained to truncate the expanded tensors. Eventually, each local tensor flows to a corresponding fixed-point tensor. One can fetch information from the fixed-point tensors, as first introduced in Ref. [25] to identify the different phases of the Ising model, where a gauge invariant quantity X is defined as 9 which represents the degeneracy of the system. It has also been used to locate the two critical points of the six-state clock model.^[18]

Taking D = 75 for an example, the calculated X is shown in Fig. 6, which equals 5 in the ordered phase and 1 in the disordered phase; while in the critical phase, it takes value in-between 5 and 1, and two abrupt jumps occur around 0.903 and 0.952, corresponding to the lower and the upper transition temperatures respectively. In this way, we can locate the transition points. Naturally, we need to perform the bond dimension scaling to acquire the converged transition temperatures, i.e., T_c1 = 0.9029(1) and T_c2 = 0.9520(1), as depicted in Fig. 7. The results agree quite well with the aforementioned MC^[15,21] and the DMRG^[20] predictions, and a brief comparison of these values is given below in Table 1 for better view.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Temperature dependence of X obtained from the fixed-point tensors, which shows two jumps around T1 = 0.903 and T2 = 0.952. The results are obtained with D = 75 as an illustration.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) The lower and upper transition temperatures with respect to bond dimensions D from X calculation. And the converged results are Tc1 = 0.9029(1), and Tc2 = 0.9520(1).

Table 1.

Table 1.

Table 1. Comparison of the critical temperatures, Tc1 and Tc2 by different methods for the five-state clock model. . Tc1 Tc2 L Tobochnik[10] (1982) 0.8 1.1 32 Borisenkoet al.[15] (2011) 0.9051(9) 0.9515(9) 1024 Kumano et al.[21] (2013) 0.908 0.944 256 Chatelain[20] (2014) 0.914(12) 0.945(17) 256 This work 0.9029(1) 0.9520(1) ∞

Table 1. Comparison of the critical temperatures, Tc1 and Tc2 by different methods for the five-state clock model. .

Briefly, we have studied the phase transitions of the ferromagnetic five-state clock model on the square lattice, using the HOTRG algorithm to investigate the thermodynamic observables with the existence of a weak external magnetic field, and the properties of the fixed-point tensor.

First, by calculating the specific heat and the correlation function, we confirm that there are indeed two phase transitions separating three phases: the low-temperature ordered phase, the high-temperature disordered phase, and the quasi-ordered phase in-between. Second, for a given bond dimension D, the upper transition temperature alone is determined by investigating the magnetic susceptibility. Through a scaling with different D’s, we obtained the converged transition temperature at T_c2 = 0.9565(7). In order to precisely determine both transition temperatures, we extracted X from the fixed-point tensors, which shows two clear abrupt jumps, corresponding to two phase transitions. Also with a bond dimension scaling, the converged transition points are located at T_c1 = 0.9029(1) and T_c2 = 0.9520(1) respectively, consistent with the MC^[15,21] and the DMRG^[20] results.